function [out,momentsModel,errorMes,model,Q] = objectFunc(paramsInputValues,setup)
% This function evaluates the objective function for GMM estimation of the RBC model
% errorMes = 0, everything is ok, 1 for reporting problems)
% Initializing the outputs
out = [];
momentsModel = [];

% For efficient starting values when computing unconditional second moments
global Var_z2old Var_z3old
if ~exist('Var_z2old','var')
    Var_z2old = [];
end
if ~exist('Var_z3old','var')
    Var_z3old = [];
end

% Check if the lower and upper bounds are violated
if any(paramsInputValues < setup.lowerBoundsValues) || any(paramsInputValues > setup.upperBoundsValues)
    errorMes = 1;
    if setup.optim == 1 || setup.optim == 3
        out = NaN;
    elseif setup.optim == 2
        out = ones(size(setup.dataMoments,1),1)*1D10;
    end
    return;
end

% We turn paramsInputValues into a struct 
for i=1:size(setup.selectParams,1)
    name = setup.selectParams(i,1);
    paramsInput.(name{1}) = paramsInputValues(i,1);
end

% Accounting for calibrated coefficients, we unold paramsInput by to get params 
params = paramsInput;
% The setdiff implies that we do not overwrite elements in params based on calibrateparams
namesCalibrate = setdiff(fieldnames(setup.calibrateParams),setup.selectParams);
for i=1:size(namesCalibrate,1)
   name = namesCalibrate(i);  
   params.(name{1}) = setup.calibrateParams.(name{1}); 
end

% Solve the model
[model,errorMes] = perturbationDSGEmodel(params,setup.orderApp,setup.selectY);
if errorMes == 1
    if setup.optim == 1 || setup.optim == 3
        out = NaN;
    elseif setup.optim == 2
        out = ones(size(setup.dataMoments,1),1)*1D10;
    end
    return;
end

% The moments in the normal distribution
ne         = size(model.eta,2);       %number of shocks
vectorMom3 = zeros(1,ne);
vectorMom4 = 3*ones(1,ne);
vectorMom5 = zeros(1,ne);
vectorMom6 = 15*ones(1,ne);

% Computing moments
sig         = 1;
maxAuto_lag = max(setup.autoLagsIdx);
mom2nd      = UnconditionalMoments_2nd_Lyap(model.gx,model.gxx,model.gss,...
              model.hx,model.hxx,model.hss,model.eta,sig,vectorMom3,vectorMom4,maxAuto_lag,Var_z2old);
Var_z2old = mom2nd.Var_z;
if setup.orderApp == 3
    mom3rd = UnconditionalMoments_3rd_Lyap(model.gx,model.gxx,model.gss,...
          model.gxxx,model.gssx,model.gsss,...
          model.hx,model.hxx,model.hss,model.hxxx,model.hssx,model.hsss,model.eta,sig,...
          vectorMom3,vectorMom4,vectorMom5,vectorMom6,mom2nd.Var_xs,mom2nd.Var_xfxf,mom2nd.Var_xs_xfxf,maxAuto_lag,Var_z3old);
      Var_z3old = mom3rd.Var_z;
end
% The moments implied by the model
if setup.orderApp == 3
    momentsModel =  momentsGMM(model,mom3rd,setup.autoLagsIdx);
else
    momentsModel =  momentsGMM(model,mom2nd,setup.autoLagsIdx);
end

% The moments implied by the model
residuals    = setup.Sw*(setup.dataMoments-momentsModel);
%Q = (dataMoments-momentsModel)'*(Sw'*Sw)*(dataMoments-momentsModel);
Q = residuals'*residuals;
if setup.optim == 1 || setup.optim == 3
    out = residuals'*residuals;
elseif setup.optim == 2
    out = residuals;
end
    

end

